Primary school temporal network

Данные

Данные взяты с этого сайта (Primary school temporal network data). Данные описывают взаимодействия учеников и учителей начальной школы в течении недели. У каждого участника эксперимента при себе было устройство, которое фиксировало каждые 20 секунд, с кем он сейчас контактирует.

Через анализ этих данных я хочу изучить на практике теорию, изложенную в книгах Network Science и Statistical Analysis of Network Data with R, и познакомиться с вариантом библиотеки igraph для языка R.

Подключаем необходимые библиотеки:

library(dplyr)
library(igraph)
library(tidyr)
library(RColorBrewer)

Читаем данные:

data1 <- read.csv("primaryschool.csv", sep="\t", header=FALSE)
colnames(data1) <- c("time", "i", "j", "classi", "classj")
data2 <- read.csv("metadata_primaryschool.txt", sep="\t", header=FALSE)
colnames(data2) <- c("i", "class", "sex")

data1 содержит информацию о том, когда (time), кто с кем (i, j) контактировал и из какого они класса или являются ли учителем (classi, classj). data2 содержит дополнительную информацию о поле учеников.

Хотим работать с взвешенным графом, поэтому посчитаем для каждой пары индивидов (ученик или учитель) число контактов (Freq) и затем удалим из данных пары индивидов, которые не контактировали друг с другом.

data1 <- data.frame(table(data1 %>% dplyr::select(i, j))) %>%
  filter(Freq > 0)

Создадим взвешенный неориентированный граф:

network <- graph_from_data_frame(data1, directed=FALSE, vertices=data2)

Проверим данные графа:

network

## IGRAPH 7522edb UN-- 242 8317 -- 
## + attr: name (v/c), class (v/c), sex (v/c), Freq (e/n)
## + edges from 7522edb (vertex names):
##  [1] 1426--1427 1426--1428 1427--1428 1426--1429 1427--1429 1428--1429
##  [7] 1426--1430 1427--1430 1428--1430 1429--1430 1426--1431 1427--1431
## [13] 1428--1431 1429--1431 1430--1431 1426--1434 1427--1434 1428--1434
## [19] 1429--1434 1430--1434 1431--1434 1426--1435 1427--1435 1428--1435
## [25] 1429--1435 1430--1435 1431--1435 1434--1435 1426--1437 1427--1437
## [31] 1428--1437 1429--1437 1430--1437 1431--1437 1434--1437 1435--1437
## [37] 1426--1439 1427--1439 1428--1439 1429--1439 1430--1439 1431--1439
## [43] 1434--1439 1435--1439 1437--1439 1426--1441 1427--1441 1428--1441
## + ... omitted several edges

Видно, что граф имеет 242 вершины и 8317 ребер.

В исходных данных нет информации о поле некоторых индивидов. Таким образом у вершин графа атрибут sex содержит значения “Unknown” заменим это значение на NA:

V(network)$sex[V(network)$sex == "Unknown"] <- NA
V(network)$sex[1:50]

##  [1] "M" "F" "M" "F" "M" "F" "F" "M" "F" "M" "M" "F" "F" "M" "F" "F" "M" "M" "F"
## [20] "F" "M" "M" "F" "M" "M" "F" "M" "F" "M" "F" "F" "M" "M" "M" "F" "F" "F" "F"
## [39] "M" "M" "F" "F" NA  "M" "F" "M" "M" "M" "M" "F"

Также вершины графа несут информацию о классе ученика или о том, что эта вершина есть учитель:

V(network)$class[1:50]

##  [1] "5B"       "5B"       "5B"       "5B"       "5B"       "5B"      
##  [7] "5B"       "5B"       "5B"       "5B"       "5B"       "5B"      
## [13] "5A"       "5A"       "5A"       "5A"       "5A"       "5A"      
## [19] "5A"       "5B"       "5A"       "5B"       "5B"       "5A"      
## [25] "5A"       "5B"       "5B"       "5A"       "5B"       "5B"      
## [31] "5A"       "5A"       "5A"       "4A"       "5A"       "5A"      
## [37] "4A"       "4A"       "5A"       "5A"       "4A"       "5A"      
## [43] "Teachers" "5A"       "4A"       "4A"       "4A"       "4A"      
## [49] "4A"       "4A"

У ребер в графе есть веса (число взаимодействий за всю неделю):

E(network)$Freq[1:50]

##  [1]  27  45   4  75 100   9  19   4   4   9  43  63  16  75  15   8  20   2  11
## [20]   4  43  12   5   4   5   7  16   3  23  44  19  62   4  42   8   6  27  13
## [39]  14  36   3  41   8  11  29   8   8   4   3   1

Анализ сети

Распределения степеней вершин

Так как у ребер есть атрибут Freq, который отражает силу связи инцидентных ребру вершин, то рассмотрим два типа распределений: распределение степеней вершин и распределение взвешенных степеней вершин (распределение сил вершин).

Добавим к вершинами атрибут degree, который будет отражать степень вершины, и wdegree, который будет отражать взвешенную степень вершины:

V(network)$degree <- degree(network)
print(paste0("the firts 10 degrees: ", toString(V(network)$degree[1:10])))

## [1] "the firts 10 degrees: 83, 47, 82, 64, 37, 71, 94, 39, 95, 76"

V(network)$wdegree <- strength(network, weights=E(network)$Freq)
print(paste0("the firts 10 wdegree: ", toString(V(network)$wdegree[1:10])))

## [1] "the firts 10 wdegree: 1109, 725, 627, 1213, 188, 1061, 1229, 508, 1910, 1000"

Посчитаем некоторые характеристики распределений степеней:

Невзвешенные степеней вершин

print(paste("Mean of sample:", mean(V(network)$degree)))

## [1] "Mean of sample: 68.7355371900826"

print(paste("Variance of sample:", var(V(network)$degree)))

## [1] "Variance of sample: 708.942217345084"

print(paste("Median of sample:", median(V(network)$degree)))

## [1] "Median of sample: 68.5"

print(paste("Min of sample:", min(V(network)$degree)))

## [1] "Min of sample: 20"

print(paste("Max of sample:", max(V(network)$degree)))

## [1] "Max of sample: 134"

Взвешенные степеней вершин

print(paste("Mean of sample:", mean(V(network)$wdegree)))

## [1] "Mean of sample: 1039.44628099174"

print(paste("Variance of sample:", var(V(network)$wdegree)))

## [1] "Variance of sample: 269076.53859607"

print(paste("Median of sample:", median(V(network)$wdegree)))

## [1] "Median of sample: 1000"

print(paste("Min of sample:", min(V(network)$wdegree)))

## [1] "Min of sample: 130"

print(paste("Max of sample:", max(V(network)$wdegree)))

## [1] "Max of sample: 2594"

Изобразим распределения степеней:

library(ggplot2)
vertices <- igraph::as_data_frame(network, what="vertices")

Невзвешенные степеней вершин

vertices %>%
  ggplot(aes(x=degree)) +
  geom_histogram(aes(y=..density..))

Взвешенные степеней вершин

vertices %>%
  ggplot(aes(x=wdegree)) +
  geom_histogram(aes(y=..density..))

Немного похожи на нормальное распределение, проверим это. Критерий Пирсона сначала делит распределение на состояния, чтобы правильно его дискретизировать. Критерий Лиллиефорса - это модификация критерия Колмогорова-Смирнова для проверки сложных гипотез. Критерий Андерона-Дарлинга - это один из критериев типа w^2. Критерий Шапира-Уилка - примерно квадрат корреляции между x и y в n.p.p.

library(nortest)

Невзвешенные степеней вершин

print(paste("Pearson", pearson.test(vertices$degree)$p))

## [1] "Pearson 0.0101508474990328"

print(paste("Lilliefors", lillie.test(vertices$degree)$p))

## [1] "Lilliefors 0.0513263641382185"

print(paste("Anderson.Darling", ad.test(vertices$degree)$p))

## [1] "Anderson.Darling 0.00278631948817735"

print(paste("Shapiro.Wilk", shapiro.test(vertices$degree)$p))

## [1] "Shapiro.Wilk 0.000596362210450844"

Взвешенные степеней вершин

print(paste("Pearson", pearson.test(vertices$wdegree)$p))

## [1] "Pearson 0.00419028224096397"

print(paste("Lilliefors", lillie.test(vertices$wdegree)$p))

## [1] "Lilliefors 0.0543648995205752"

print(paste("Anderson.Darling", ad.test(vertices$wdegree)$p))

## [1] "Anderson.Darling 0.001974323336201"

print(paste("Shapiro.Wilk", shapiro.test(vertices$wdegree)$p))

## [1] "Shapiro.Wilk 0.000114865431336385"

Критерий Лилифорса с уровнем значимости 0.5 не отверг гипотезу о нормальности для обоих распределений, остальные критерии отвергли.

Далее будем работать с взвешенными степенями, так как это более естественно для взвешенного графа.

Корреляция степеней вершин

Посчитаем среднюю взвешенную степень для каждой вершины: (anwd = average neighbor weighted degree) и изобразим результат на графике:

V(network)$anwd <- knn(network, weights=E(network)$Freq)[[1]]
vertices$anwd <- V(network)$anwd

vertices %>%
  ggplot(aes(x=degree, y=anwd)) +
  geom_point() +
  geom_smooth(method="lm") +
  xlab("Weighted Vertex Degree") +
  ylab("Average Neighbor Weighted Degree")

Можно видеть положительную корреляцию. Интересно посмотреть на этот же график, но с разделением по классам. Снизу слева видна отдельная группа людей, интересно, не учителя ли это.

colors <- c("deeppink1","cyan4",
            "brown1","olivedrab4",
            "cyan1","deeppink4",
            "olivedrab1","royalblue3",
            "gold1","gold4",
            "black")

vertices %>%
  ggplot(aes(x=degree, y=anwd, color=class)) +
  geom_point() +
  geom_smooth(method="lm", se=FALSE) +
  xlab("Weighted Vertex Degree") +
  ylab("Average Neighbor Weighted Degree") +
  scale_color_manual(values=colors)

Нет, это не были учителя.

Посчитаем без разделения по классам корреляцию Пирсона для степеней вершин, то есть наклон регрессионной прямой, изображенной на предыдущем графике, а также корреляцию Спирмена:

cor(vertices$degree, vertices$anwd, method="pearson")

## [1] 0.5823455

cor(vertices$degree, vertices$anwd, method="spearman")

## [1] 0.5750678

Больше 0.5. Сеть можно назвать ассортативной, то есть вершины в ней предпочитают связываться с узлами, похожими на них самих. В нашем случае похожесть отражается во взвешенных степенях вершин.

Центральность вершин

Найдем топ 10 центральных вершин, основываясь на степенях. Я предполагаю, что самыми социальными являются учителя. Помимо взвешенной степени будем использовать степень близости, степень посредничества и степень влиятельности.

V(network)$closeness <- closeness(network, weights=1/E(network)$Freq)
vertices$closeness <- V(network)$closeness
V(network)$betweenness <- betweenness(network, weights=1/E(network)$Freq)
vertices$betweenness <- V(network)$betweenness
V(network)$evcent <- evcent(network, weights=1/E(network)$Freq)$vector
vertices$evcent <- V(network)$evcent

Взвешенная степень

(vertices %>% arrange(-wdegree))[1:10,]

##      name class sex degree wdegree     anwd  closeness betweenness    evcent
## 1695 1695    1B   F    103    2594 88.66577 0.10481716        2233 0.7349942
## 1697 1697    1B   M    116    2532 88.89297 0.10949262        3363 0.8352363
## 1665 1665    1B   M    123    2447 82.07233 0.10895842        1983 0.8854687
## 1698 1698    1B   M     99    2442 85.66298 0.11131071        3995 0.6211609
## 1613 1613    2A   M    111    2249 93.96398 0.10784996        4263 0.7322550
## 1666 1666    1B   M    110    2234 93.70054 0.10638859        1557 0.8345163
## 1912 1912    1B   M    111    2227 93.40323 0.10922895        2947 0.8182644
## 1675 1675    1B   F    104    2197 82.02822 0.10573275        2158 0.7990344
## 1558 1558    3B   M    101    2119 84.52336 0.08618615        2023 0.6358605
## 1890 1890    2B   M    120    2096 75.57061 0.09525315        1309 0.8259771

Степень близости

(vertices %>% arrange(-closeness))[1:10,]

##      name class sex degree wdegree     anwd closeness betweenness    evcent
## 1698 1698    1B   M     99    2442 85.66298 0.1113107        3995 0.6211609
## 1697 1697    1B   M    116    2532 88.89297 0.1094926        3363 0.8352363
## 1912 1912    1B   M    111    2227 93.40323 0.1092290        2947 0.8182644
## 1665 1665    1B   M    123    2447 82.07233 0.1089584        1983 0.8854687
## 1663 1663    1B   F     55    2060 92.16214 0.1085594        1037 0.3535092
## 1613 1613    2A   M    111    2249 93.96398 0.1078500        4263 0.7322550
## 1688 1688    1B   M     86    1893 88.88906 0.1074014        2665 0.6063156
## 1666 1666    1B   M    110    2234 93.70054 0.1063886        1557 0.8345163
## 1675 1675    1B   F    104    2197 82.02822 0.1057327        2158 0.7990344
## 1708 1708    2B   M    118    1839 82.06743 0.1052624        2652 0.8983541

Степень посредничества

(vertices %>% arrange(-betweenness))[1:10,]

##      name class sex degree wdegree     anwd  closeness betweenness    evcent
## 1613 1613    2A   M    111    2249 93.96398 0.10784996        4263 0.7322550
## 1698 1698    1B   M     99    2442 85.66298 0.11131071        3995 0.6211609
## 1697 1697    1B   M    116    2532 88.89297 0.10949262        3363 0.8352363
## 1912 1912    1B   M    111    2227 93.40323 0.10922895        2947 0.8182644
## 1664 1664    1B   M     85    1503 87.64205 0.10356201        2925 0.6376254
## 1837 1837    4A   M     96    1582 71.72819 0.09651677        2862 0.6037813
## 1688 1688    1B   M     86    1893 88.88906 0.10740143        2665 0.6063156
## 1708 1708    2B   M    118    1839 82.06743 0.10526243        2652 0.8983541
## 1533 1533    4A   M     96    1498 75.77637 0.09481543        2319 0.6269478
## 1695 1695    1B   F    103    2594 88.66577 0.10481716        2233 0.7349942

Степень влиятельности

(vertices %>% arrange(-evcent))[1:10,]

##      name class  sex degree wdegree     anwd  closeness betweenness    evcent
## 1551 1551    3B    M    134    1514 88.12219 0.08072743         343 1.0000000
## 1761 1761    1A    M    128    1197 84.16291 0.08792606         410 0.9985940
## 1780 1780    3A    M    129    1534 89.13233 0.08526683         946 0.9392089
## 1552 1552    3B    F    123    1684 84.65499 0.08290923         671 0.9245904
## 1673 1673    1B    M    124    1568 80.01658 0.09590599           3 0.9043510
## 1708 1708    2B    M    118    1839 82.06743 0.10526243        2652 0.8983541
## 1665 1665    1B    M    123    2447 82.07233 0.10895842        1983 0.8854687
## 1765 1765    1B    F    118    1845 77.50081 0.09871025         585 0.8844436
## 1730 1730    4A <NA>    108    1192 73.14933 0.07999364         909 0.8821280
## 1560 1560    3B    F    100    1434 86.36820 0.08013365         475 0.8811915

Нет, оказалось, что учителя не являются самыми социальными. Интересно, что среди самых социальных индивидов много парней. Посмотрим на примере класса 4B и его учителя, насколько он менее социален детей из класса. Изобразим подграф в форме целевой диаграммы (чем ближе точка к центру тем больше значение параметра индивида). Так как имеем дело с плотным графом, ребра много информации нести не будут. Окрасим мальчиков и девочек в синий и розовый цвета соответственно, а учителя отметим черным.

subgraph <- induced_subgraph(network, V(network)[class=="4B" | name=="1521"])
E(subgraph)$inv.Freq <- 1/E(subgraph)$Freq
V(subgraph)[class=="4B" & sex=="M"]$color <- "dodgerblue"
V(subgraph)[class=="4B" & sex=="F"]$color <- "pink"
V(subgraph)[class=="Teachers"]$color <- "black"
A <- as_adjacency_matrix(subgraph, sparse=FALSE)
library(network)
g <- network::as.network.matrix(A)
library(sna)

Взвешенная степень

sna::gplot.target(g, V(subgraph)$wdegree,
   circ.lab = FALSE,
   circ.col="skyblue", usearrows = FALSE,
   vertex.col=V(subgraph)$color,
   edge.col="darkgray")

Cтепень близости

sna::gplot.target(g, V(subgraph)$closeness,
   circ.lab = FALSE,
   circ.col="skyblue", usearrows = FALSE,
   vertex.col=V(subgraph)$color,
   edge.col="darkgray")

Cтепень посредничества

sna::gplot.target(g, V(subgraph)$betweenness,
   circ.lab = FALSE,
   circ.col="skyblue", usearrows = FALSE,
   vertex.col=V(subgraph)$color,
   edge.col="darkgray")

Cтепень влиятельности

sna::gplot.target(g, V(subgraph)$evcent,
   circ.lab = FALSE,
   circ.col="skyblue", usearrows = FALSE,
   vertex.col=V(subgraph)$color,
   edge.col="darkgray")

На графиках можно видеть, что учитель - самый несоциальный индивид, так как расположен почти всегда дальше всех от центра. Можно долго исследовать характеристики центральности и по-разному визуализировать результаты, но остановимся на этом.

Плотность сети

Вычислим плотность сети:

edge_density(network)

## [1] 0.2852097

Вычислим плотности подсетей для каждого класса, они по отдельности должны быть плотнее всей сети школы:

class <- unique(V(network)$class)

for(i in class) {
  subgraph <- induced_subgraph(network, V(network)[class == i])
  print(paste("density of", i, ":", edge_density(subgraph)))
}

## [1] "density of 5B : 0.978260869565217"
## [1] "density of 5A : 0.974025974025974"
## [1] "density of 4A : 0.995238095238095"
## [1] "density of Teachers : 0.644444444444444"
## [1] "density of 3B : 0.991341991341991"
## [1] "density of 4B : 0.944664031620553"
## [1] "density of 2A : 0.972332015810277"
## [1] "density of 1B : 0.996666666666667"
## [1] "density of 2B : 0.969230769230769"
## [1] "density of 1A : 0.984189723320158"
## [1] "density of 3A : 0.99604743083004"

Так и выходит, логично, что дети внутри класса контактируют друг с другом больше.

Изобразим матрицу сопряженности для всей сети, чем темнее точка тем выше степень взаимодействия двух индивидов:

edges <- as_long_data_frame(network) %>% 
  mutate(from = as.character(from), to = as.character(to))

tmp <- data.frame(to = edges$from, from = edges$to, Freq = edges$Freq, 
                 to_class = edges$from_class, from_class = edges$to_class)

tmp <- rbind(edges %>% select(from, to, Freq, to_class, from_class), tmp)

tmp %>%
  ggplot(aes(from, to)) +
  geom_raster(aes(fill=Freq)) +
  scale_fill_gradient2(low="gray100", mid="gray50", high="gray0") +
  theme(legend.position="none") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        panel.background = element_blank())

Матрица сопряженности демонстрирует плотность сети, отсутствие точки означает, что два индивида не взаимодействовали. Изобразим матрицу сопряженности для одного класса, например, 4B:

tmp %>% filter(from_class == "4B", to_class == "4B") %>%
  ggplot(aes(from, to)) +
  geom_raster(aes(fill=Freq)) +
  scale_fill_gradient2(low="gray100", mid="gray50", high="gray0") +
  theme(legend.position="none") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        panel.background = element_blank())

Можно видеть, что матрица сопряженности намного плотнее для конкретного класса.

Длины путей

Посчитаем среднее расстояние между индивидами в сети.

mean_distance(network)

## [1] 1.732451

В принципе ожидать, что в школе, где есть учителя, через которые связанны все классы, будут длинные пути нелогично. Найдем максимальную длину:

max(distances(network))

## [1] 3

Коэффициент кластеризации

Посчитаем коэффициенты кластеризации для всей сети (вероятность того, что смежные для случайной вершины вершины будут связаны):

transitivity(network)

## [1] 0.4797899

Посчитаем коэффициенты кластеризации для каждой вершины:

V(network)$transitivity <- transitivity(network, type="local")
V(network)$transitivity[1:50]

##  [1] 0.5301205 0.5735430 0.5401987 0.5188492 0.6846847 0.5098592 0.4468085
##  [8] 0.7287449 0.4759239 0.5382456 0.8487903 0.4723684 0.4992743 0.5184552
## [15] 0.5369674 0.4525799 0.4247843 0.5079365 0.4796099 0.5519126 0.4661503
## [22] 0.6023222 0.5536232 0.4285714 0.5412262 0.4519597 0.5459877 0.5782313
## [29] 0.5480769 0.4957895 0.5245372 0.4541204 0.4386343 0.4667269 0.5256140
## [36] 0.5361492 0.4359514 0.5662879 0.7300945 0.4629870 0.5268293 0.5113459
## [43] 0.4692308 0.4929577 0.9249012 0.5078631 0.5102041 0.4875776 0.4798246
## [50] 0.8000000

vertices$transitivity <- V(network)$transitivity

Проверим топ 10 индивидов с наибольшим коэффициентом кластеризации:

(vertices %>% arrange(-transitivity))[1:10,]

##      name class sex degree wdegree     anwd  closeness betweenness     evcent
## 1750 1750    5B   M     21     255 68.07451 0.05415851           0 0.06326984
## 1715 1715    2B   F     22     334 50.17066 0.06586774           0 0.09601607
## 1524 1524    4A   F     23     412 73.09466 0.07232686           0 0.09579966
## 1441 1441    5B   M     32     270 69.17407 0.05499871           0 0.19418689
## 1538 1538    4A   F     26     429 66.94406 0.06842667           0 0.06773123
## 1897 1897    2B   F     30     502 63.01793 0.07351401           0 0.07640525
## 1803 1803    4B   M     24     299 71.23746 0.07218493           0 0.11684756
## 1609 1609    2A   F     20     388 39.06959 0.06426053           0 0.14595304
## 1661 1661    1B   F     33    1301 75.52037 0.09263242         282 0.13438010
## 1737 1737    3A   M     47     729 88.07133 0.07635733          31 0.28066311
##      transitivity
## 1750    1.0000000
## 1715    0.9913420
## 1524    0.9249012
## 1441    0.8487903
## 1538    0.8000000
## 1897    0.8000000
## 1803    0.7789855
## 1609    0.7631579
## 1661    0.7537879
## 1737    0.7335800

Это ученики с небольшими степенями.

Анализ соседей учителей

Зная, что в школах России у начальных классов есть по одному главному учителю, который преподает почти все дисциплины, интересно проверить, будем ли мы наблюдать ту же картину в иностранной школе. Посчитаем соседей для каждого учителя в каждом классе:

for(i in V(network)[class == "Teachers"]) {
  print(paste0("neighbors of the teacher with id ", V(network)[i]$name, ":"))
  print(table(neighbors(network, i)$class))
}

## [1] "neighbors of the teacher with id 1521:"
## 
##       1B       3A       3B       4A       4B       5A       5B Teachers 
##        1        1        1        4       23        1        2        7 
## [1] "neighbors of the teacher with id 1650:"
## 
##       1B       2A       2B       5B Teachers 
##        6       12       26        5        5 
## [1] "neighbors of the teacher with id 1653:"
## 
##       1A       1B       2A       3A       4A       4B       5A       5B 
##        1        6        1        1       21        3        3        3 
## Teachers 
##        7 
## [1] "neighbors of the teacher with id 1668:"
## 
##       1A       1B       2A       3A       3B       4A       4B       5A 
##        3        3        1        1        5        1        1       22 
##       5B Teachers 
##        7        5 
## [1] "neighbors of the teacher with id 1709:"
## 
##       1B       2A       2B       3A       3B       4A       4B       5A 
##        3        1        1       18       22        1        2        2 
##       5B Teachers 
##        3        6 
## [1] "neighbors of the teacher with id 1745:"
## 
##       1A       1B       2A       2B       3A       3B       4A       5A 
##        9       25        3        2        3        1        1        2 
##       5B Teachers 
##        3        7 
## [1] "neighbors of the teacher with id 1746:"
## 
##       1A       3A       3B       4A       4B       5B Teachers 
##        1       23       15        3        2        2        5 
## [1] "neighbors of the teacher with id 1753:"
## 
##       1A       1B       2A       4A       4B       5B Teachers 
##       20        5        2        3        2        2        4 
## [1] "neighbors of the teacher with id 1824:"
## 
##       3B       4A       5A       5B Teachers 
##        1        1        5       24        5 
## [1] "neighbors of the teacher with id 1852:"
## 
##       1A       1B       2A       2B       3A       3B       4A       4B 
##        3        7       23       10        5        1        1        3 
##       5B Teachers 
##        1        7

Видно, что чаше у одного учителя подавляюще много контактов с учениками одного класса, но также встречаются учителя, которые активно контактируют с учениками двух классов.

Визуализация сети

Протестируем разные способы визуализации сети.

Методы визуализации

Начнем с круговой визуализации. Так как в сети много вершин, то визуализация всех данных через круг приведет к тому, что на ней будет ничего не видно. Но визуализировать отдельный класс не имеет смысла, так как получится почти полный граф. Посмотрим как выглядят связи для топ 30 учеников с наименьшими степенями свободы.

set.seed(0)

igraph_options(vertex.size=8, 
               vertex.label.cex=0.5, 
               vertex.label.dist=1.5, 
               vertex.label.degree=3*pi/2,
               edge.width=1)
subnetwork <- subgraph(network, (vertices %>% filter(class != "Teachers") %>%
                                   arrange(wdegree))[1:30,]$name)
plot(subnetwork, layout=layout_in_circle)

Продолжим работать с той же выборкой учеников. Попробуем разные методы визуализации, позже визуализируем всю сеть. Сеть:

plot(subnetwork, layout=layout_on_grid)

Сфера:

plot(subnetwork, layout=layout_on_sphere)

Случайное расположение вершин:

plot(subnetwork, layout=layout_randomly)

Алгоритм Дэвидсона и Харела:

plot(subnetwork, layout=layout_with_dh)

Алгоритм Фрухтермана и Рейнгольда:

plot(subnetwork, layout=layout_with_fr)

Сигма алгоритм:

plot(subnetwork, layout=layout_with_sugiyama)

В общем, в igraph есть много разных методов отрисовки графов, не будем все тестировать.

Декорирование графа

Перейдем к визуализации всей сети, только теперь будет декорировать граф, используя возможности библиотеки igraph. Каждый класс будем окрашивать в свой цвет. Вершины учеников изобразим как круги, а вершины учителей - как квадраты. Размер вершины будет пропорционален ее степени. Ширина ребра пропорционально ее весу.

V(network)$label <- V(network)$name

V(network)$shape <- "circle"
V(network)[class=="Teachers"]$shape <- "square"

V(network)[class == "1A"]$color <- "deeppink1"
V(network)[class == "1B"]$color <- "cyan4"
V(network)[class == "2A"]$color <- "brown1"
V(network)[class == "2B"]$color <- "olivedrab4"
V(network)[class == "3A"]$color <- "cyan1"
V(network)[class == "3B"]$color <- "deeppink4"
V(network)[class == "4A"]$color <- "olivedrab1"
V(network)[class == "4B"]$color <- "royalblue3"
V(network)[class == "5A"]$color <- "gold1"
V(network)[class == "5B"]$color <- "gold4"
V(network)[class == "Teachers"]$color <- "black"

E(network)$color <- "black"

set.seed(1)
plot(network, 
     vertex.size=sqrt(V(network)$wdegree/50), 
     edge.width=E(network)$Freq/400,
     vertex.label=NA,
     layout=layout_with_fr)

На этом графике хорошо заметно как учителя занимаются с одним или двумя классами, не более.

Граф более высокого уровня

Теперь объединим вершины учеников по классам и построим граф для более высокого уровня. Каждая вершина соответствует своему классу, размер вершины отображает число учеников в нем. Толщина ребра отражает степень взаимодействия между классами. В центре расположим учителей, а классы расположим по кругу.

classes <- as.factor(V(network)$class)
network.c <- contract(network, classes, 
                      vertex.attr.comb=list(name="ignor",
                                           class="first",
                                           sex="ignor",
                                           degree="ignor",
                                           transitivity="ignor",
                                           label="ignor",
                                           shape="first",
                                           color="first"))
V(network.c)$label <- V(network.c)$class
V(network.c)$class.size <- as.vector(table(V(network)$class))
network.c <- simplify(network.c,
                      edge.attr.comb=list(Freq="sum",
                                          color="ignor"))
E(network.c)$color <- "azure4"

set.seed(1)
l <- layout_as_star(network.c, center = 11)
plot(network.c,
     vertex.size=V(network.c)$class.size,
     edge.width=E(network.c)$Freq/500,
     vertex.label.dist=3,
     vertex.label.color="black",
     layout=l)

График хорошо показывает, ученики каких классов активно друг с другом взаимодействуют. Если сравнить этот график с тем, что выше, можно увидеть, что толщина ребер графа для классов хорошо отражает пересечение (близость друг к другу) классов в графе. Вероятно, у этих классов есть общие предметы или кружки. Видно, что толстые ребра между классами сверстников.

Эгоцентрическая сеть

Построим эгоцентрическую сеть, в центре которой, например, будет учитель с наименьшим индексом. В такой сети отражены соседи выбранной вершины и все связи между ними. Учителя с индексом 1521 поместим в центр и изменим его цвет на белый.

network.e <- make_ego_graph(network, order=1, nodes=V(network)[name=="1521"])
V(network.e[[1]])$label <- ""
V(network.e[[1]])[name=="1521"]$color <- "white"

set.seed(0)
plot(network.e[[1]], 
     vertex.size=sqrt(V(network.e[[1]])$wdegree/6), 
     edge.width=E(network.e[[1]])$Freq/60,
     layout=layout_with_fr)

Видно, что две наибольшие группы, с которыми учитель 1521 взаимодействует, - это другие учителя и класс 4B.

Резюме

Книга Network Science представляет из себя что-то среднее между научно-популярной литературой и плохим учебником. Книга Statistical Analysis of Network Data with R великолепна, концентрация интересной и полезной информации с приведенным кодом для более легкой реализации поражает. Я ранее пользовался библиотекой igraph, когда программировал на Python, использовать ее в R оказалось приятнее, скорее всего, потому что я люблю R… В Rmd файле уже больше полтысячи строк кода и комментариев, и он долго исполняется, так что на этом анализ завершаю. При этом осталось очень много нетронутых мной базовых и более продвинутых тем в области анализа сетевых данных!